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I.NTKOL  uC'i  luN 

In  recent  years  there  has  been  a vide  acceptance  of  linear 
tractuio  meclaric^  insulting  in  tin  develi.pi  cut  o!  now  structural 
alloys  having  I ibli  1 1 ,i<  tare  toughm  and  u.aint airing  yield  strength 
close  to  pro'  1 u-  w \t  1 - . 

However,  plastic  it>  plays  a major  role  in  the  application  of 

these  material*  cithti  in  thin  cross  sections  or  under  mixed  mode 

eonditiors.  A In.  in  some  cases,  tr  meet  the  ASTI*  requirement  for 

plane  strain  fracture  toughness  testing,  the  speamns  required  are 

too  large  fei  iceror.iodl  testing.  To  alleviate  some  of  these  problems 

a number  of  retie. ds  have  been  proposed,  c.g.,  Irwin's  equivalent  'F.lastic 

Crack  l.engtl  Kell':  (iack  Opening  Displacement",  Rice's  Path 

Independent  .1- Irie^ral'  «<ru  Non- linear  htiergy  Methods  proposed  by 

fc1 

I lebnuitz  aiid  his  ioW'  tiers  , the  last  two  being  quite  promising, 
lienee  it  is  r.roe  - .ary  tr  model  the  plastic  condition  near  the  crack  tip 
a.  accurately  j i 'Eil.lt  . 


1 Irwin,  ti.  R.,  Fracture  Testing  of  lii gli- Strength  Sheet  Materials  Under 
Conditions  Aj ; noprinte  lor  Stress  Analysis,  Naval  Research  Laboratory, 
Ppt.  548b.  July  l'.'fO. 

■-Kells,  A.  A.,  hn«-taHe  Crack  Priq  ; gat  i en  ii  Meta  1 , Proc . Conf.  Crack 
Propagation,  (inrlicld,  hngland,  19h0,  p.  120. 

‘Rico,  .1.  i .,  A r t V.  Independent  Integial  and.  the  .Appi  oximate  Analysis 
of  ‘train  Core  t t t..t  ’•  n h>  Notches  and  Cracks,  Trais.  Am.  Soc . Mech. 
Ingr'-.,  ..'oiin.il  of  Applied  Mechanics,  1 r 6 c , p.  A " 9 . 

'*!  1 1 i - , d . , , 1 . I , and  Iiebowitz,  I'.,  On  True  ture  Toughness  in 

the  Nnnlir.tat  i n.gi  , Ingr.  f-rnct  . Mech.,  Vul.  ",  11)75,  p.  491. 
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V 
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1 


Lx;. 


Jk 


X- 


In  this  paper,  we  implement  higher  order  isoparametric  elements 
(quadratic  and  cubic)  in  NASTRAN's  piecewise  linear  (plasticity) 
module.  By  judicious  choice  of  intermediate  grid  points,  and  using 
proper  constraints,  we  develop  elastic  and  elastic-plastic  singular 
elements. 

Specifically,  the  elastic  singular  cubic  element  embodying  the 
square  root  (l/Zr)  singularity  is  constructed  by  placing  the  midside 
nodes,  adjacent  to  the  crack  tip,  at  l/9th  and  4/9th  locations.  The 
plastic  singular  element  is  constructed  for  the  Ramberg-Osgood  type  of 
material  with  zero  hardening  exponent  (ideally  plastic  material)  using 
the  'Sliding  Node  Concept'  of  Barsoum.^ 

'Sliding  Nodes'  are  simply  achieved  by  collapsing  one  side  of  an 
element  and  surrounding  the  crack  tip  with  these  elements,  so  that  the 
crack  tip  has  multiple  independent  nodes  at  one  physical  location 
which  slide  with  respect  to  each  other  during  deformation,  due  to 
loading.  The  proper  order  for  plastic  singularity  (i.e.,  1/r)  is 
achieved  by  locating  the  adjacent  midside  nodes  at  l/9th  and  4/'9th 
of  the  length  of  the  side  of  the  element,  as  done  for  the  elastic 
element . 

After  a brief  review  of  the  theory  proving  the  existence  of  crack 
tip  singularities,  we  discuss  the  implementation  of  these  elements  in 
NASTRAN  as  user  dummy  elements.  The  results  of  the  analysis  are 
compared  to  a Prandtl  slip-line  field  solution. 

-’Barsoum,  R.  S.  , Triangular  Quarter-Point  Elements  as  Elastic  and 
Perfectly-Plastic  Crack  Tip  Element,  Int.  J.  Num.  Meth.  Engrg.  , Vol. 
1977,  p.  85. 
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Many  general  purpose  finite  element  codes  as  well  as  advanced 
versions  of  NASTRAN  may  have  these  elements.  Hence,  the  method  may  be 
i,m te  accessible  to  many  users. 

CRACK  TIP  SlNCULARITlhS 

Consider  the  path  integral  J developed  by  Rice"  ’ , 


J = J (Wdy  - T * |H.  ds) 


(1) 


where  l\ 
vector s 
a crack 


is  the  strain  energy  density,  r and  u traction  and  displacement 
on  the  path  F.  Using  a circular  path  of  radius  r surrounding 
tip  (1)  reduces  to. 


J 


if 

r / {W  cosG  - x 
-IT 


3u 

37 


■}d0 


(2) 


The  terms  in  {...}  ir.  above  are  of  the  form: 

(stress)  (strain) , 

hence  for  the  nonvanishing  contribution  to  J (which  is  identical  to 
energy  release  rate  for  the  elastic  case),  we  have 


o.  . £.  . = 0 (— ) a s r -*•  0 . 
1J  ij  tj 


(3) 


liquation  (3)  is  quite  familiar  for  the  elastic  case  for  which  stress  and 
'train  each  have  a singularity  of  the  order  of  one  half  at  the  crack  tip. 


’Rico,  J.  R.,  A Path  Independent  Integral  and  the  Approximate  Analysis 
c.f  ft  rain  Concentration  by  Notches  and  Cracks,  Trans.  Am.  Soc . Mech. 
Ingrs.,  Journal  of  Applied  Mechanics,  1968,  p.  379. 

^Rice,  .J.  R.  and  Rosengren,  G.  F.,  Plane  Strain  Deformation  Near  a 
Crack  Tip  in  a Pcwer-I.aw  Hardening  Material,  J.  Mech.  Phys.  Solids, 
Vol.  16,  p.  1. 


Now  consider  the  Ramberg-Osgood  type  of  material  given  by 


T 


T = CY  = Z~  T > Tor  y < yq 
'0 

(4) 

n 

T = Tq^7-)  » for  Y 

' n 

(5) 

where  x = /T72S~~S~  , y = , and  xQ,  yQ 

strain  in  shear  and  n is  the  hardening  exponent. 

(3),  outside  the  elastic  range,  we  have 

-n 

n+1 

= 0(r  ) 

J 1_ 

J *n 


are  yield  stress  and 
From  (4) , (5)  and 


(6) 


From  (6)  we  have  the  familiar  elastic  case  for  n = 1.  However  when 
n = 0,  which  is  the  case  of  ideally  plastic  material,  we  have  from  (6) 


°ij  = 0(r°) 


(7) 


Ei;j  = 0(r_1) 

indicating  a singularity  of  order  one  for  the  strains. 

The  existence  of  such  singularities  for  quadratic  elements  have 
been  given  in  ref.  5,7  . In  the  next  section  we  briefly  outline  the 
case  of  the  cubic  elements. 


^Barsoum,  R.  S. , Triangular  Quarter-Point  Elements  as  Elastic  and 
Perfectly-Plastic  Crack  Tip  Element,  Int . J.  Num.  Meth.  Engrg. , 

Vol . 11,  1977,  p.  85. 

Hussain,  M.  A.,  Lorensen,  IV.  E.,  and  Pflegl , G.,  The  Quarter-Point 
Quadratic  Isoparametric  Element  as  a Singular  Element  for  Crack 
Problems,  NASTRAN  Users'  Experiences,  NASA  TM-X-3428,  Oct.  1976,  p.  419. 


SINGULARITIES  OF  CUBIC  ELEMENTS 


Following  the  notation  of  ref.  8 , the  geometry  of  a 12-point  cubic 
element  is  mapped  into  a normalized  square  in  (£,n)  plane  (-1  < £ < 1 , 

-1  <_  ri  < 1)  through  the  transformation 


12 

x = l N^CC.rOxi  , 
i = l 

12 

y = I Nj(C,n)yi  , 

i = l 


(8) 


where  the  shape  function  is  given  by 


Ni  = 2k  (1  + cci)(1  + r'rii)['10  + + n2) ] [-10  ♦ 9(5?  + n2)] 

o 1 

♦ — (i  + 55. Hi  + 9nn.)(i  - n2)(i  - n2) 

256  1 1 1 

♦ ^ (i  + nn.)d  + 9«J(1  - C2)  (1  - 5?)  , (9) 


(x^,  y^  and  5^,  are  the  grid  points.) 

Collapsing  the  quadrilateral  element  as  shown  in  Figure  1 and 
placing  the  midside  nodes  at  l/9th  and  4/9th  location,  we  have 

X1  = X10  = X11  = x12  = °*  x2  = x9  = hi/9’  X3  = X8  = 4h/9  * 


X 


4 


x 


5 


x 


6 


h * ho  ’ hi  ' h2  ‘ °-  h ’ ->'9  ‘ -*/9'  h * ' -u/9 

U = ~y~  = >'s  s "y6  = ‘ Vi  ' K 

o , 

“Zienkiewicz , C.O.,  The  Finite  Element  Method  in  Engineering  Science, 
McGraw  Hill,  London,  1971. 


Substituting  (10)  into  equation  (8)  we  have 

x - - (1  + O2 
4 

y = ^ (i  + 5)2  dl 

4 

1/2 

Any  point  at  a radial  distance,  r = (x2  + v2)  ' , from  the  crack  tip 

is  given  by 

r = | (l  + 02[(|)2  + n2]]/z 


(1+0  = vr 


{JU£)2  + h2]1/2l1/2 


The  Jacobian  [J]  is  given  by 


[J]  = 


3x  3y 

W W 

3x  9y 
3i7  cirj- 


jd  + o yici  + o 


id  + o: 

4 


and  the  determinant  is 


det  | J | = ^ (1  + O3 


For  the  inverse  functions,  we  have 


K 

9]1 

dX 

3x 

3n 

9y 

9y 

2 -4n 

h Cl  + O h Cl  + O' 


F(i  + O' 


The  displacement  components  of  the  point  (5,n)  for  an  isoparametric  trans- 


formation are, 


12 


u = l Ni(C,n)ui 


i=l 

12 


(16) 


v=  l n. (5,n)v. 
i=l 


The  derivatives  of  u,  v 

with 

respect 

to  5,  n 

are 

3u 

12 

Y 

9Ni 

3u 

12 
_ Y 

85 

= L 

i=l 

-5TUi 

* 

- L 
i=l 

8n 

3v 

85 

12 

8N. 

8v 

’ 8n 

12 

3N. 

= l 

i=l 

— — v. 
85  1 

II 

H- 

II 

l 

8n 

ui 


(17) 


where 


3N 

ar-  = 2k  (1  + nni)["10  + 9^i  + n2i)l(-io^i  + 95^  + 185  ♦ 275^) 

81 


35  256 


+ TT-  (!  - nf)5id  + 9nn.)(i  - n2) 

♦ d + nni)  (l  - 52)(95t  - 25  - 275^*) 


(18) 


9N4  , , 

— = — (.1  + 55.)  [-10  + 9(52  + n2)](-i0n.  + 9n.52  ♦ I8n  ♦ 27n.n2) 

dr)  256  1 1111  1 

O 1 

+ 256  (1  + ^iU1  ' ni)(9T1i  ' 2n  " 

+ — (l  - 52)n.  (l  + 955. )d  - 52) 

256  1 1 1 

9 

Substituting  for  nodal  values  and  collecting  terms,  using  MACSYMA, 
equations  (17)  become 


9MACSYMA:  Math  Lab  Group,  MIT  Laboratory  for  Computer  Science  (Symbolic 
Manipulation  System),  November  1975. 


f =a0+al(l+5)  +a2(l+0 


^ = b0  + b2  (1  + 5)  ♦ b2(l  + O2  * b3(l  ♦ £)3  (19) 

[9(-3u12  + 3uu  - u1Q  ♦ u7  - 3u6  ♦ 3u5  - u4  ♦ u^n* 

9(U12  + U11  ' U10  + U7  " U6  - U5  + U4  ' Ul)n2 

+ (27u12  ' 27U11  - 35U10  + % ' 9U8  ' U7 

+ 27u,  - 27u_  + u + 9u  - 9u  + 35u  )n 

6 5 4 3 2 1 

('9U12  " 9U11  ’ 35U10  + 9U9  - 9U8  • U7  + 9U6 

+ 9Ug  - - 9uj  + 9u0  - 35Uj)]  (20) 

- 5ug  + 4ug  - u?)(l  + n)  - (u4  - 4u3  + 5u2  - 211^(1  - r,)  ] 


1Q  - 3ug  + 3ug  - u?)(l  + n)  - (u4  - 3u3  + 3u2  - i^)  (1  - n)  ] 


7(3u12  • 3ull  * “10  - “l’1'  ' 1S(U12  * U11  - U10  - Ul)n 


* U1  - 27U12  * 27“ll  - U10] 


(21) 


[27 (3Uj2  “ 3u^  ^10  ~ u7  ” 3U5  + u4  •*  11^)0 


18(u12  + un  - u1Q  + u7  - u6  - u5  + u4  - Uj)n2 


+ (*27u12  + 27U11  + 35U10  ‘ 72U9  + 36U8  + U7 


- 27u,  + 27u  • u,  - 36u  + 72u  - 35u  )] 


(22) 


p 

I 


b9  = — (2\Jjq  - 5uq  + 4ufi  - u7  + u;1  - 4u,  + 5u?  - 2u,) 


16 


= - — (ujq  - 3uq  + 3ug  - Uy  + U4  - 3ug  + 3 u}  - Uj) 


The  derivatives  3v/9£,  9v/9n  are  the  same  except  for  replacing  by  v^. 
The  derivatives  of  u with  respect  to  x,y  are  obtained  from 


9u  _ 9u  95  9u  9r[ 


9x  9£  3x  3p  3x 


4t|bo 


2ap  - 4pbi  1 


h(l  + 5)2  h(l  + 5) 


+ ^ (2aL  - 4nb2)  + r-  (1  + 5)  (2a-,  - 4pb3) 


(23) 


9u  9u  9£  9u  9p 


3y  9^  9y  9n  9y 


4b0  4b 1 4b2  4b3 

+ + * — (1  ♦ f) 

«.(!  + 5) 2 m + o 1 i 


(24) 


Similar  expressions  are  obtained  for  9v/9x  and  3v/3y  with  u replaced 

1 


by  v in  n's  aid  L's. 
j 


The  stresses  and  strains  are  singular  when  the  Jacobian  determinant 
vanishes  at  £ = -1.  From  (23),  (24)  and  (12),  the  singularity  is  0(l/r) 
if  b(J  jt  0 and  is  0(l//r)  if  bQ  = 0.  A careful  study  of  (21)  indicates  that 
b()  depends  on  the  displacements  of  nodal  points  at  the  crack  tip.  If 
the  nodal  points  at  the  crack  tip  are  tied  together,  i.e., 


U1  = U10  = U11  = U12 


and 


V1  = V10  = V11  = V12 


(25) 


then  hn  = 0 and  the  strain  field  has  the  inverse  square  root  of  r 
singularity,  the  correct  singularity  of  linear  fracture  mechanics.  On 
the  other  hand  if  the  nodal  points  at  the  crack  tip  are  allowed  to  move 
independently  to  one  another,  the  strain  field  has  the  (1/r)  singularity, 
a characteristic  of  perfect  plasticity. 

0 


_ J 


NASTRAN  IMPLEMENTATION 


The  NASTRAN  implementation  for  the  quadratic  element  follows  the 
steps  outlined  in  section  6.8.3.12  of  reference  12.  The  following 
routines  require  modification:  PLA1 , which  creates  the  ECPT's  and 
EST's  for  the  linear  and  non-linear  elements;  PLA31  and  PLA32,  which 
recover  stresses  for  the  non-linear  elements;  and  PLAYBD,  PLA41  and 
PLA42  which  control  generation  of  the  updated  stiffness  matrix.  The 
following  new  routines  are  required:  PSDUM1,  a driver  for  stress  data 
recovery  in  PLA3;  PSDM11  and  PSDM12,  phase  I and  II  stress  recovery 
routines;  PKDUM1 , a driver  for  stiffness  generation  for  the  non-linear 
elements;  PKDM11  and  PKDM12,  stress  recovery  routines  which  generate 
stresses  for  the  computation  of  the  non-linear  material  matrix;  and 
PKDM1S,  the  stiffness  matrix  generation  routine  for  non-linear  elements. 
The  two  driver  routines,  PSDUM1  and  PKDUM1  can  be  modelled  after  the 
corresponding  routines  for  the  QUAD1  element.  The  remaining  routines 
are  modifications  of  the  stiffness  and  stress  recovery  routines  required 
for  rigid  format  1,  statics.  The  major  modifications  to  switch  from 
statics  to  piecewise  linear  include  changing  the  labelled  common  areas, 
building  the  non-linear  material  matrix^0  and  calculating  incremental 
stress  rather  than  total  stresses. 

Hussain,  M.  A.,  Lorensen,  IV.  E.,  and  Pflegl,  G.,  The  Quarter-Point  Quad- 
ratic Isoparametric  Element  as  a Singular  Element  for  Crack  Problems, 
NASTRAN  Users'  Experiences,  NASA  TM-X-3428,  Oct.  1976,  p.  419. 

10The  NASTRAN  Theoretical  Manual,  Editor,  MacNeal,  R.  H.,  NASA  SP-221, 

Sep.  1970,  p.  104. 

12The  NASTRAN  Programmer's  Manual,  NASA  SP-223(01),  Sep.  1972. 
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NUMERICAL  EXAMPLE 


Consider  the  problem  of  small  scale  yielding.  The  problem  is 
governed  by  the  elastic  field  at  points  far  away  from  the  crack  tip  and 
asymptotically  has  the  elastic  singular  field.  Near  the  crack  tip  we 
have  the  plastic  zone.  This  is  schematically  represented  in  Iigurc  2.  The 
plane  strain  slip  line  field  is  also  shown. 

The  problem  is  modelled  in  a fashion  similar  to  Barsoum's^.  The 
geometry  is  shown  in  Figure  3.  The  crack  tip  elements,  1-12,  are  the 
singular  elements  which  can  either  be  quadratic  or  cubic  elements.  For 
the  symmetric  case  the  corner  nodes  of  the  elements  are  placed  on 
concentric  semi-circles,  0 <_  0 <_  it,  at  tt/12  intervals,  of  radii,  r = 0, 

.5,  1.0,  1.625,  1.52,  22,  2.52,  32 , 42,  S.52. 

The  method  of  solution,  for  the  plastic  problem,  is  based  on 
Swedlow's  piece-wise  linear  analysis  and  is  well  documented  in  the 
NASTRAN  theoretical  manual. ^ 

The  procedure  for  the  present  problem  is  accomplished  via  two 
rigid  formats.  The  static  rigid  format  is  first  used  to  obtain  the 
stress  distribution  and  the  equivalent  stresses  at  the  integration 
point  (£  = n = 0)  for  the  elastic  increment.  This  solution  is  performed 
with  all  the  collapsed  nodes  at  the  crack  tip  having  the  same  displacement 

^Barsoum,  R.  S.,  Triangular  Quarter-Point  Elements  as  Elastic  and 
Perfect ly-Plastic  Crack  Tip  Element,  Int.  J.  Num.  Meth.  Engrg., 

Vol . 11,  1977,  p.  85. 

lf)The  NASTRAN  Theoretical  Manual,  Editor,  MacNeal,  R.  H.,  NASA  SP-221, 

Sept.  1970,  p.  104. 
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vector  (see  equat  ion  25).  This  is  accompli  shed  with  multipoint 
constraints.  The  outermost  nodes  are  subjected  to  the  displacements 


governed  by  Westergaard ' s solution,  with  K = 1, 


" -4cf)1/2  (f^-cos  0, 

v 2G  2tt  1 + v 


cos  0/2 
sin  0/2 


where  E = 30  x 10  psi  and  v = .3.  The  value  of  2K^  is  established 
from  the  elastic  solution  based  on  the  yield  stress  (aQ)  of  20  x 10‘^ 
psi  for  the  highest  stressed  element.  For  the  plastic  analysis  the 
stress-strain  curve  is  provided  with  the  above  constants  and  yield  strain 
at  .2%  and  hardening  exponent  n * .3  (this  should  be  close  to  zero  for 
perfect  plasticity).  The  nodes  at  the  crack  tip  are  then  released  for 
sliding  in  order  to  obtain  1/r  singularity  at  the  crack  tip.  The  load 
is  incremented  by  K^/4  till  the  plastic  zone  has  reached  the  first  layer 
of  elements. 

Preliminary  results  of  the  problem  are  indicated  in  Figures  4a-c 
and  compared  with  those  of  ref.  11  . From  the  static  solution  it  was 
found  that  the  inception  of  yielding  occurs  at  0 = 68°  compared  to  the 
theoretical  value  of  0 = 70° . 


In  Figures  4a-c  we  have  also  plotted  the  slip  line  (plane  strain) 
solution  for  comparison.  The  plastic  zone  also  corresponds  well  with 
ref.  11. 


Hutchinson,  .1.  W. , Singular  Behavior  at  the  End  of  a Tensile  Crack 
in  a Hardening  Material,  .1.  Mech.  Phys.  Solids,  1968,  Vol.  16,  p.  13. 
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CONCLUSION 

Higher  order  isoparametric  elements  can  be  effectively  used  for 
modelling  singular  elastic  as  well  as  plastic  problems  that  arise  in 
the  field  of  fracture  mechanics.  The  procedure  in  obtaining  these  do 
not  require  any  special  crack  tip  elements  but  are  simply  constructed 
by  adjusting  the  adjacent  nodes  at  proper  locations  and  proper 
constraints.  The  locations  of  these  nodes  should  be  adhered  to  as 
closely  as  possible  for  stable  answers.  Since  many  general  purpose 
finite  elements  may  have  these  elements  in  their  library  the  method, 
for  crack  problems,  may  be  accessible  to  many  users. 
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COLLAPSED  TRIANGULAR  ELEMENTS 
AROUND  A CRACK  TIP 


FIGURE  3.  FINITE  ELEMENT  IDEALIZATION  OF  THE  CRACK  TIP  NEAR  FIELD. 


FIGURE  4.  (a)  SHEARING  STRESS  DISTRIBUTION  FOR  THE  SINGULAR  ELEMENTS  COMPARED  TO  REF.  (11)  AND  A SLIP-LINE 
SOLUTION  (PRELIMINARY  RESULTS). 
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